Draft version November 16, 2011 

Preprint typeset using 1^1^^ style emulatoapj v. 02/07/07 



ORIGIN OF MULTIPLE NUCLEI IN ULTRALUMINOUS INFRARED GALAXIES 

HiDENORi Matsui\ Takayuki R. Saitoh^"', Junichiro Makino^'^'^, Keiichi Wada^, Kohji Tomisaka^-^, Eiichiro 
KoKUBO^'^, HiROSHi Daisaka"\ Takashi Okamoto", & Naoki Yoshida' 

(Received- Accepted) 
Draft version November 16, 2011 



ABSTRACT 

Ultraluminous infrared galaxies (ULIRGs) with multiple (> 3) nuclei are frequently observed. It has 
been suggested that these nuclei are produced by multiple major mergers of galaxies. The expected 
rate of such mergers is, however, too low to reproduce the observed number of ULIRGs with multiple 
nuclei. We have performed high-resolution simulations of the merging of two gas-rich disk galaxies. 
We found that extremely massive and compact star clusters form from the strongly disturbed gas disks 
after the first or second encounter between the galaxies. The mass of such clusters reaches ~ 10* M©, 
and their half-mass radii are 20 — 30 pc. Since these clusters consist of young stars, they appear to 
be several bright cores in the galactic central region (~ kpc). The peak luminosity of these clusters 
reaches ~ 10 % of the total luminosity of the merging galaxy. These massive and compact clusters 
are consistent with the characteristics of the observed multiple nuclei in ULIRGs. Multiple mergers 
are not necessary to explain multiple nuclei in ULIRGs. 

Subject headings: galaxies: starburst — galaxies: interactions — galaxies: evolution — method: 
numerical 



1. INTRODUCTION 

Ultraluminous and luminous infrared galaxies 
(ULIRGs/LIRGs) have strong infrared (IR) luminosities 
of LiR > 10 ^^ Lp, and 10" L© < Ltr < lO^^ L©, 
respectively (jSanders fc Mirabel 1 19961 ) . Their intense 
IR luminosities are mainly due to starburst activities, 
and the contribution of active galactic nuclei (AGN) to 
IR luminosity seems to be rather l imite d , i.e., ^ 20 % 
(iRisaliti et all 120061: iFarrah et ahl 120071: iNardini et al] 
I2008D. The starburst activities in (U)LIRGs are likely 
to be triggered by mergings of galaxies, in particu- 
lar gas-rich galaxies, since observations suggest that 
(U)LIRGs have complex and disturbed morphologies 
fSander s fc~Mirabellll996[ ). In the /-band (F814W filter) 
images of 100 sampled (U)LIRGs at 0.05 < z < 0.20 
taken with Hubble Space Telescope (H ST), about 20 % of 
ULIRGs have m ultiple (> 3) nuclei (jBorne et al.l 120001 : 
ICui et al.l[200l . The fraction increases to more than 
80 % , if all probable galaxies are included (jBorne et al.l 
[2OOOI) . 

The origin of the multiple nuclei has been thought 
to be multiple major merger event s . In order to ex- 
plain multiple mergers, iBorne et all ()2000[ ) argued that 
the progenitors of (U)LIRGs, were compact groups of 
galaxies. However, the evolution timescale of compact 
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groups is very long Hq^) (jAthanassoula et aI1ll997f l. 
Therefore, the probability that two galaxies merge in 
a compact group and the merging galaxy has double 
cores is Ppair ^ t^mg/H^^ ~ 100 Myr/10 Gyr - 0.01, 
where imcrg is the merging timescale of the galactic cores. 
The probability that another galaxy merges with such 
merging galaxy in the compact group and the merging 
galaxy has triple cores is Pmuiti = Ppair x j^j^:^ x ppair ~ 

3 X 10~^, where Ng (~ 5) is the typical number of galax- 
ies in a compact group (|Mendes de Oliveira &: HicksonI 
Il991h . Then, the number density of ULIRGs with mul- 
tiple cores is 0g x Pmuiti ~ 3 x 10~^ Mpc^'^, where 0g 
(~ 10"'* Mpc~'^) is the number density of galaixes in 
a compact group (iMendes de Oliveira &: HicksonI 119911 : 
iRibeiro et all I1994D . Since the number density of 
ULIRGs is 10-^ Mpc"^ (jSanders et al.l[2003h . the frac- 
tion of multiple merger to ULIRGs is (3 x 10~^)/10~^ ^ 
0.03. Thus, it is unliky that multiple merger is a major 
cause of the observed multiple nuclei in (U)LIRGS. An 
alternative explanation for the origin of multiple nuclei 
in (U)LIRGs seems necessary. 
Pre vious numer ic al sim u lations of merging galax- 



ies (iMihos et al.l [l99l iMihos fc Hernauist 
i Barnes fc Hernquis tI Il996t iKazantzidis et al. 
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Cox et al.. 2 006: Di Matteo et al. 2007; Narav anan et al.l 
2009 HCovington et al.l l20in[r have shown that starbursts 
occurred during the merging process. However, these 
simulations failed to reproduce formation of star clusters 
as observed in many in teractin g and me rging galaxies 
(jWhitmore fc Schweize? 1995; M engel et a l. 2008), since 
their spatial and mass resolutions {e.g., 100 pc and 
^ 10^ -^0) were not sufficiently high to distinguish star 
forming region and formation of star clusters whose sizes 
are less than 100 pc. In addition, previous simulations 
did not allow the interstellar medium (ISM) to cool 
below ~ 10* K. This is another reason why formation 
of star clusters was not reproduced in the previous 
simulations. 
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In order to reproduce formation of star clust ers, high 
resolu tion simulations have been attempted. iLi et alJ 
()2004|) have performed the simulation of a merging 
galaxy using Tree+SPH code GADGET. In this sim- 
ulation, mass and spatial resolutions were 10 pc and 
6.6 X 10^ -Mq, respectively. They assumed an isothermal 
ISM and used sink particles, that absorb their surround- 
ing gas, to represent clusters. They have shown that a 
number of massive clusters form in the merging process. 
In this simulati on, the most ma s sive c luster has mass 
of 7.8 X 10"^ Mq. iBournaud eFaLl (|2008[ ) have performed 
the simulations of merging galaxies using sticky particles, 
that collide inelastically, instead of solving hydrodynam- 
ics. In these simulations, mass and spatial resolutions 
were 32 pc and 7 x 10^ M©, respectively. They have also 
shown formations of massive star clusters with masses of 
10^"^ Mq. These simulations are, however, unrealistic 
in a sense that dynamical evolution of star clusters them- 
selves cannot be properly followed due to t he limitation 
of the sticky and sink particle methods. ISaitoh et al.l 
(|2009l) improved the spatial and mass resolutions (5- 
20 pc and 10^"'' Mq) and the ISM model that is allowed 
to cool to 10 K. These simulations showed that the be- 
havior of the multiphase ISM in the merging galaxies is 
considerably altered and the formation of shock-induced 
star clusters is naturally reproduced (Saitoh et al. 2009, 
1201 It ). There are several high resolution simulations 
of merging galaxies, which resolve the low temperature 
gas (300- 500 K) using adaptive mesh refinemen t (AMR) 
methods (jKim et al.l 120091 : iTevssier et al.ll2010D . These 
simulations also showed the difference in the behavior of 
the multiphase I SM from that of ISM used in the previ- 
ous simulations. iKim et aD ()2009l ) considered the merger 
of low-mass galaxies (~ 1.8 x 10"^0 Mq). Several spiky 
peaks of star formation rate were seen in thier simulation 
although the most prominet starburst was found at the 
beginning the simualtions and the significant fraction of 
the ISM was blown away by energetic wind before the 
merging event. Thus, it would be hard to investigate 
detailed evolu t ion of the ISM in merging galaxies. In 
ITevssier et al.l ()2010l ). the polytropic equation of state 
was used instead of solvi ng the energy equat ion, which is 
essentially different from lSaitoh et all (|2009l ) and present 
paper. The validity of this approximation is unclear for 
understanding the evolution of the ISM in the merging 
galaxies. 

High-resolution images of the local (U)LIRGs ob- 
tained by the integral fi eld spectroscopy usi ng Willia m 
Herschel Telescov e (iGarcia- Man'n et all I20q9alf 



and VLT-VIMOS ('Alonso-Herrero et a l.l 120091 1201 
iMonreal-Ibero et al . 2010; Rodriguez Zau rin et al1l2010( ) 
showed that (U)LIRGs generally have very complex 
structures, such as Ha bright knots, rings, and tidal 
tails. As in the case of formation of star clusters, the 
reason why these structures have not been reproduced 
in numerical simulations might be just the inadequate 
treatment of ISM and limited resolution. Thus, high 
resolution simulations resolving multiphase ISM are 
essential to comprehend the complex structures in 
(U)LIRGs. 

We have performed high-resolution simulations of 
merging galaxies with sufficiently high spatial resolution 
and cooling function of ISM that covers a wide range of 
temperature (10 K < T < lO'* K). In this paper, we fo- 



cus on the origin of the observed ULIRGs with multiple 
nuclei. The detailed evolution of merging galaxies will 
be described in forthcoming papers. 

The structure of this paper is as follows. We first de- 
scribe the numerical methods and models in f}2] Numeri- 
cal results are shown in SjH In 21 'we compare our results 
with observations. Conclusions and discussion are pre- 
sented in fjs] 

2. METHODS 
2.1. Simulation Setup 

The model paramete rs of the initia l disk galaxy are the 
same as those used in ISaitoh et al.l ([2009) . The masses 
of the dark matter halo and old stellar and gas disks 
are 1.1 x 10" Mq, 5.1 x 10^ Mq, and 1.2 x 10^ Mq, 
respectively. The gas fraction in the disk is ^ 20% of 
the total disk mass. The gas fraction is set to be slightly 
higher than that of local spiral galaxies since the gas frac- 
tion decreases due to star formations during the isolated 
phase for 1 Gyr before interactions (see below). We also 
assume the halo gas component of which total mass is 
1.1 X 10^ Mq. It has the same distribution as the dark 
matter halo. The initial gas temperatures are 10^ K in 
the disk and 10^ K in the halo. The scale radii of the 
stellar and gas disks are 4 kpc and 8 kpc, respectively. 
The dark-matter halo and the stellar disk are expressed 
by A^-body particles, whereas the gas disk is expressed 
by smoot hed particle hydro dynamics (SPH) particles. 

Unlike ISaitoh et aD ()2009f ). before starting the simula- 
tions, we let an isolated disk galaxy evolve for 1 Gyr 
in order to stabilize the gas component. After this 
evolution, the amount of disk and halo gas in each 
galaxy is 1.8 x 10^ Mq. The SFR during this period is 
^ 0.5 Mq yr~^. During the first 1 Gyr, the disk galaxy 
holds a quasi-steady state and is free from a global in- 
stability due to gravity. The evolution of the isolated 
disk is different from that of a gas rich disk at high- 
z. Simulations of the gas-rich disk, in which gas mass 
fraction of a disk is > 50 %, have shown that massive 
clum ps form from gravitational instabill i ty of the gas rich 
disk jNoguchil 119991: llmmeh et ahl 12004 IBournaud eFall 
|2007() . In contrast to the simulations of rnerging of such 
clumpy disks at high-z bv IBournaud et all (|2011[ ). we fo- 
cus on the merging process of two local disk galaxies. 

We let the pre-evolved two galaxies collide in a 
parabolic orbit with the pericentric distance of 7.5 kpc. 
The simulations start with the initial separation of 
75 kpc, where the separation distance is measured be- 
tween the mass centers of the two galaxies. 

We have performed nine runs. Three of them are 
high-resolution runs (i?xx in Tab. [1]) and the rests are 
low-resolultion runs (Lxx)- The initial masses of SPH, 
star, and dark matter particles are the same. They are 
7.5 X 10^ Mq and 3 x 10'' M©, for Hxx and Lxx runs, 
respectively. Subscripts P, T, and R denote prograde, 
tilted, and retrograde spin axes of two galaxies. We per- 
formed six possible combinations of the spin orientations 
in low resolution and only PP and TT cases in high res- 
olution. The gravitational softening length, e, is set to 
be 20 pc for eight runs except for -ffxT.Spc run in which 
the softening length is 5 pc. 

The orientation of the disk axis is specified by two an- 
gles, i and uj. Here, i is the inclination and uj is the 
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argument of pericenter. They are defined as follows. We 
use the coordinate system in which the orbital plane of 
two galaxies is in the x-y plane. The Laplace-Runge-Lenz 
vector of the orbit is = (1,0,0). The spin axis of a 
galaxy is then given by eg — (sin a; sin i, cos a; sin z, cos i). 
These angles are defined relative to the orbit of each 
galaxy. Thus, in the case of "TT" runs, spin axes of 
two galaxies are parallel to each other. 

2.2. Numerical Method 

Numerical simulations were performed by an A^- 
body/SPH code ASURA (Sai toh et all l2009l) that can op- 
tionally utilize GRAPE ( Sugimoto et all |1990D in or- 
der to accelerate calculation of gravitational force. We 
adopted a time-step limiter (jSaitoh fc Ma kino 2009) for 
the time-integration of SPH particles with individual 
time-steps, which keeps the differences in time-steps of 
neighboring particles to be small enough to handle the 
strong shocks correctly. We also used the FAST scheme 
(jSaitoh fc Makinall2 010') that allows self-gravitating SPH 
particles to use different time-steps for integrations of 
hydrodynamics and gravity. This method accelerates 
simulations without losing the accuracy of the time- 
integration. The leapfrog method is used for the time 
integrations of both gravity and hydrodynamics. These 
improvements allow us to follow the formation of cold 
and dense gas clumps and their expansion by supernovae 
(SNe) feedback without numerical problems. 

For the ra diative cooling and heating, we follow the 
treatment of ISaitoh et al.l (|2008[ ). in which the gas is al- 
lowed to cool down to 10 K. The treatment of star for- 
mation and SN fee dback are also the same as those in 
ISaitoh et all ()2009[ ). When an SPH particle satisfied (1) 
high number density (tih > 100 cm^"^), (2) a low temper- 
ature [T < 100 K), and (3) a collapsing region (V-t> < 0), 
the SPH particle spawns a star particle of which mass is 
one-third of the original SPH particle mass. We assume 
simple stellar polulation approximation for newly formed 
star particles. Salpeter's initial m ass function wit h mass 
range of 0.1 - 100 Mq is adopted ()Salpeteij|19"55l ). Stars 
with mass > 8 Mq in each star particle explode as Type 
H SNe and release 10^^ erg of thermal energy per one SN 
into the surrounding ISM. 

3. MULTIPLE NUCLEI IN ULIRGS 

Initially, the two galaxies move on the parabolic orbits 
and approach to each other. At around t = 450 Myr, 
they reach the pericenter (first encounter). Since or- 
bital angular momenta of their main bodies are con- 
verted into the disk internal spin, they los e their orbital 
angular moment a (Barnes 1992; iHernqu ist 1992, 1993; 
iMihos fc Hernquisti 119961) and their trajectories deviate 
from the parabolic orbits. The main bodies approach 
to each other again (second encounter). Finally, after 
several encounters, their cores completely merge. 

Figures [T] and [5] show the evolution of distribution of 
gas and stellar particles for runs i?pp and Htt- Gas 
disks are strongly disturbed at the first encounter (the 
panels at t — 430, 463, and 596 Myr) and at the second 
encounter (the panels at t = 830 and 863 Myr) in run 
Hpp and at the second encounter (the panels at t = 863 
and 896 Myr) in run Htt- Because of the strong per- 
turbation, gas is compressed by large-scale shocks. The 
dense gas radiatively cools and becomes gravitationally 



unstable. As a result, m assive star clusters form (see also 
ISaitoh et al. l [201(il[20Tlh . 

Figure |3| shows the synthesized /-band images from 
our merger simulations after the massive star clusters 
formed. Here, the /-band luminosity emitted from newly 
formed stars was calculated using population synthesis 
code PEGASE ()Fioc fc Rocca-Volmerangelll999f) . and ef- 
fects of dust absorption and re-emission were neglected. 
The two top panels show the high resolution runs, and 
the six lower panels show the low resolution runs. In all 
images, multiple compact sources appear in the central 
regions of galaxies. Two of them with arrows are the pro- 
genitor galactic cores, and the others are newly-formed 
massive star clusters. The absolute magnitudes of these 
clusters in /-band are Mj < —17.0, which are compara- 
ble to those of the galactic central cores. The effect of 
dust extinction is discussed in §4.2. 

Since the cores form through the gravitational instabil- 
ity of strongly perturbed gas disks, they do not contain 
dark matter and consist only of young stars. On the other 
hand, the galactic luminous nuclei contain dark matter 
and a large amount of old stars that have formed before 
the merger event. The properties of newly formed cores 
are, therefore, different from those of the galactic nuclei. 
In addition, such cores are more massive than usual star 
clusters with 10®~^ Mq, by an order of magnitude. We, 
therefore, call the newly formed cores "hypermassive star 
clusters" . Since the hypermassive star clusters form in 
the galactic central region kpc), they are different 
from tidal dwarves that are formed in tidal arm s (e.g., 
Barnes fc Hcrnquist 1992; Bournaud et al.ll2008D . 

The reason why they are extremely massive can be 
explained as follows. The total mass of the gas which 
becomes gravitationally unstable is very large (~ 1.0 x 
10^ Mq) because of the strong large-scale gravitational 
and liydrodynamical disturbances from another galaxy at 
the encounter phase. The strong disturbances deform the 
gravitational potential into strongly non-axisymmentric 
shape. Therefore, the rotation of the disk becomes 
ineffective in stabilizing long-wavelength perturbations. 
Even though the gas disk fragments to a number of small 
clumps with ^ 10^~^Mq, these small clumps are still 
gravitationally bound to each other and eventually merge 
to form larger clusters. 

Figure HI shows the snapshots during the formation of 
the most massive cluster from 900 Myr to 913 Myr in run 
Htt- Multiple gas clumps formed from the disturbed 
gas disk merge with each other and become one massive 
cluster. In other words, the merging of smaller clumps 
within the short timescale causes quick growth of the 
mass of star clusters. As a result, hypermassive star 
clusters form. 

In Fig. |5l we show the evolution of the masses of these 
hypermassive star clusters (upper panels) , and their dis- 
tances from the center of the galaxy (lower panels). Here, 
the definition of the distance is the smaller of the two dis- 
tances from the two galactic centers. The clusters first 
grow quickly, reaching about a half of the final mass in 
around 20 Myr. After that, the growth of the mass slows 
down, and the clusters fall to the center of the galaxy 
due to the dynamical friction. As a result, the merger 
remnant has a very compact and luminous core similar 
to those observed in a large number of elliptical galaxies 
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TABLE 1 
The model parameters. 



Run 


1 L oj 


uj^ [dee;! 


[dee;] 


uj^ [deffl 

2 L 6J 


lit. |^J.^J.(^J 




star 




e' [dcI 














7.5 X 10^ 


514476 


1841042 


27720000 


20 




-109 


-30 


71 


-30 


7.5 X 10^ 


514476 


1841042 


27720000 


20 


Lpp 












3.0 X 10* 


133696 


447994 


6930000 


20 









180 




3.0 X 10-* 


133696 


447994 


6930000 


20 


Lbs. 


180 




180 




3.0 X lO'* 


133696 


447994 


6930000 


20 


Lpx 







71 


30 


3.0 X 10-* 


133696 


447994 


6930000 


20 




-109 


30 


180 




3.0 X 10-* 


133696 


447994 


6930000 


20 


Z/XT 


-109 


-30 


71 


-30 


3.0 X 10* 


133696 


447994 


6930000 


20 


-f^TT.Spc 


-109 


-30 


71 


-30 


7.5 X 10^ 


514476 


1841042 


27720000 
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^ The disk inclination for the galaxy 1. ^ The argument of the pericenter for the galaxy 1. The disk inclination for the galaxy 2. The 
argument of the pericenter for the galaxy 2. The mass of SPH, star, and dark matter particles. ^ The number of SPH particles. ^ The 
number of star particles. ^ The number of dark matter particles. ' The gravitational softening length. 



(e.g., iKormendv et al.ll2009D . 

Figure [S] shows the evolution of the SFR and bolo- 
metric luminosity coming from the stars in the merging 
event. The effects of dust absorption and re-emission are 
not taken into account. If we take them into account, 
the IR luminosity would become higher since ultraviolet 
and optical photons are absorbed by dust and re-emitted 
in the IR regime. Therefore, we expect that the bolo- 
metric luminosity of our runs is a good indicator of the 
IR luminosity. In all runs, the maximum SFR reaches 
20 M0yr~'^, and the luminosity reaches 10^^ Lq, which is 
about an order of magnitude higher than that before the 
encounter. The figure shows that the evolution of bolo- 
metric luminosity reflects the SFR. The peak luminos- 
ity is comparable to those of LIRGs (jSanders fc Mirabel 
Il996l ). If the progenitor galaxies are more massive and/or 
the initial gas fraction is much higher than those used in 
this simulation, the peak of SFR would become higher 
and the merger galaxy would become ULIRGs with in- 
frared emission > lO^^L©. AGNs would also provide 
an additional IR luminosity, although our simulations do 
not take into account the effect of AGNs. The influence 
of AGNs to the total IR luminosity is typically ~ 20 % 
of the total one (iRisali ti et all [20061: iFarrah et "all 120071: 
iNardini et al.ll2008[) . The peak mass-to-light ratio of our 
run, M/ L, is less than 0.1 and is comparable to that of 
ULIRGs (iGolina et all [20051 : iHinz fc Riekg [20061) . The 
period of the high luminosity phase is ~ 200 Myr. The 
hypermassive star clusters after their formation contain 
a large amount of dusty gas and young stars as shown 
in Fig. m Therefore, a strong infrared emission from the 
hypermassive star clusters is likely to be observable. 

In order to study the dependence on the numerical res- 
olution of the formation of hypermassive star clusters, 
we have performed simulations with different mass reso- 
lutions. In both high and low resolution runs, hypermas- 
sive star clusters with ~ 10^ Mq are formed as shown 
in Fig. [3l Their formation process is the same as that of 
the high resolution runs. In addition, the time evolution 
of SFR is also similar to that of high resolution runs. In 
both runs, the peak of SFR is 20 - 30 Mq yr~\ and 
the duration time of the active star formation is about 
200 Myr. 

4. COMPARISON WITH OBSERVATIONS 

4.1. Surface Brightness and Spatial Distribution of 
Multiple Nuclei 



Here, we compare our num erical resul t s with observa- 
tions of ULIRGs with nuclei. ICui et ahl (|2001l ) observed 
nuclei with absolute magnitude from —17 to —21 in I- 
band. The /-band luminosity of the identified putative 
nuclei is typically around a fe w percent of far -infrared 
luminosity of their host galaxy ()Gui et al.ll2001t ). 

We estimated the absolute /-band magnitude of hy- 
permassive star clusters formed in our simulations using 
PEGASE. Some hypermassive star clusters have luminosity 
comparable to the observed putative nuclei. In the upper 
panels of Fig. [31 there is one hypermassive star cluster 
with AIj ~ —17.4 in run //pp, and there are three hyper- 
massive star clusters with A// ~ —17.7, —17.1, and —17.1 
in run /^tt, respectively. These /-band absolute magni- 
tudes are comparable to those of the galactic cores, which 
are —18.5 and —17.5 in run Hpp and —17.7 and —15.7 
in run //tt- These results indicate that the galactic nu- 
clei and hypermassive star clusters are not distinguish- 
able from their luminosities. The /-band (bolometric) 
luminosity of hypermassive star clusters is a few percent 
10%) of the bolometric luminosity of the host galaxy 
with ^ 10^^ Lq. These characteristics are in good agree- 
ment with the observations. 

We also compare the spatial distribution of the ob- 
served putative nuclei with that of simulated hyper- 
massive star clusters. If we select relatively compact 
ULIRG s, the av erage separation between nuclei is about 
1 kpc (|Cui et al . 2001). In our simulations, the sepa- 
ration between hypermassive star clusters and galactic 
cores depends on the evolution phase. However, the typ- 
ical separation is a few kpc in the most luminous phase 
of the galaxies. Thus, the spatial distribution also agrees 
with the observations. 



4.2. Dust Extinction 

In this subsection, we investigate the effect of the dust 
extinction in /-band. Here, we use the high spatial reso- 
lution run, //tt,5pc7 in which the softening length is 5 pc 
for mock obscrvatations, since the detailed structure of 
the ISM is important to estimate the extinction. In this 
run, similarly to run //tt, hypermassive star clusters 
with > 10* Mq form after the second encounter. The 
masses of the first, second, third most massive star clus- 
ters are ~ 1.6 x 10® Mq, 1.1 x 10* A/©, and 6.1 x 10^ Mq, 
respectively. 

The flux, taking dust extinction into account, from 





10-^1 10^10=^ [Mg pc"^] 1 lO' 10^ [Mq pc"^] 10-'' 1 lO' [Mg pc"^] 

Fig. 1. — Snapshots of run Hpp. The upper color, middle gray-scale, and bottom color panels show surface density of gas, old stars, and 
newly formed stars that are born after t = 0, respectively. In each upper panel, the number in the right-bottom corner displays simulation 
time (the unit is Myr). The white line in the right-bottom corner of each panel shows length of 5 kpc. 



a star is estimated by F = Fq exp(— r^). Here, 

i 

Fq is the flux expected without dust extinction, 
and Ti is optical depth for e ach SPH particl e and 
n = o-H.T A^H.i {Zi/Zo,) (|Bekki &: Shioval I2OOOI : 
iKobavashi et al.ir2010[ ). where ^d, ctr,!, -^H,i, and Zi rep- 
resent the dumpiness parameter, the dust cross-section 
per H atom at /-band for the Milkyway dust model, col- 
umn number density of atomic H of a SPH particle, and 
the metallicity of a SPH particle, respectively. The SPH 
particles between the star and the observer are used for 



calculating the optical depth, r. The dumpiness param- 
eter is introduced in order to represent the dumpiness of 
the dust distribution in the ISM. It is assumed to be 0.15, 
0.5, and 1.0. The value of (u = 0.15 is used in Lya emit- 
ters (jKobavashi et al.ll2010l ). whereas that of = 1-0 is 
the case without the dust dumpiness under the spatial 
resolution. Here, we do not estimate the "dumpiness" 
of SPH particles along the line of sight directly, since 
the "real" dumpiness comes from much small er scale 
in th e ISM. We adopt ch,! = 4.0 x 10"^^ cm^ (|Drainel 
|2003( ). The column density of an SPH particle is given by 
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10-^1 10^10=^ [Mg pc"^] 1 lO' 10^ [Mq pc"^] 10-'' 1 lO' [Mg pc"^] 

Fig. 2. — Same as Fig. [T]but run Htt- The snapshots are projected to the disk plane. 



-^H,i = 'm'iW{r^ h), where rrii and W{r^h) are the mass 
and 2-dimensional kernel function of the SPH particle, 
respectively, where r is the 2-diniensional distance be- 
tween the SPH and star particles in the projected plane, 
and h is the kernel size of the SPH particle. 

Figure [7] illustrates the effect of dust extinctions. The 
upper row shows the time evolution of the /-band map 
without extinction. The middle row represents the ex- 
pected observation taking dust extinction with = 0.5 
into account for run i?TT,5pc- Here, the map is projected 
to the orbital plane of the galaxies. After formation of 
hypermassive star clusters, strong dust extinctions take 



place in the clusters since there is a large amount of 
dusty gas. The dust extinction reduces the flux from 
the formed hypermassive star clusters. The extinction 
of the clusters in /-band is A{I) > 5 magnitude in the 
qd — 0.5 case as shown in the bottom panels of Fig. [71 
In the models of qd = 0.15 and qd ~ 1.0, the extinc- 
tions are A{I) > 4 and A{I) > 6, respectively. The 
strong extinction continues till 970 Myr (the left and the 
middle panels), which corresponds to the epoch about 
50 Myr after the formation of the clusters. In this dusty 
phase, multiple core structures are buried by dust. After 
990 Myr (the right panels), the clusters become optically 
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Fig. 3. — Simulated HST WFPC2 7-band (F814W filter) images from our simulations. The color represents magnitude per square 
arcsecs. The left and right panels in the top row show the results of runs -ffpp and i^xTi respectively. These two runs are simulated with 
high resolution. Six panels in the bottom two rows show the results of low resolution runs. We ignore the effect of dust extinction (see 
§4.2.). 



thin because of consumption of dusty gas by star for- 
mation and escape of dusty gas from the cluster by SN 
feedback. As a result, the multiple cores appear. The 
extinction becomes Aj < 1 magnitude. On the other 
hand, the galactic central region is still buried by dust 
due to the continuous gas accretion. In this phase, the 
extinction does not depend strongly on the dumpiness 
parameter since dusty gas is restricted to the central re- 
gion of the cluster rather than distributed in the broad 
region of the cluster. 

The timescale in which hypermassive star clusters be- 
come optically thin in /-band is less than the ULIRG 
lifetime (~ 200 Myr as shown in Fig. |6|). Therefore, it is 
possible for the clusters to be observed as ULIRGs. Note, 
however, that the estimate of the dust extinction depends 
strongly on the distribution of the dust and that resolu- 
tions of current simulations are insufficient to resolve the 
"real" structure of the ISM. Futhcr high resolution sim- 



ulations are necessary to involve the dumpiness of the 
dust distribution in the ISM. 

4.3. The Fraction of ULIRGs with Multiple Nuclei 

The fraction of ULIRGs with multiple nuclei strongly 
depends on the wavelength used for observation. In the 
an alysis of HST J-ban d dat a of ^ 100 ULIRGs samples 
by iBorne et all ()2000[ ) and ICui et all (|2001l) . the frac- 
tion of t he ULIRGs with multiple nuclei is ~ 20%. In 
additi on. ICui et all (l2001l) have analyzed the nine sam- 
ples of ISurace et"aLl (| 19981) in /-band. Two ULIRGs are 
found to have multiple nuclei, i.e., the fra ctio n is consis- 
tent w ith the results of iBorne et al.l (|2000( ) and lCui et al.l 
(|200l . 

On the other h and, iBushouse et al.l (|2002l ) and 
IVeilleux et al.l ()2002l ) have claimed that the fraction of 
ULIRGs with multiple nuclei is less than 5%, based on 
their analysis of the data o bserved in near I R ba nds. 
Using the HST //-band data, IBushouse et al.l (|2002D re- 
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Fig. 4. — The formation of the most massive star cluster in run i/xT (the black curve in the right panels of Fig. [5j. The upper panels 
show the evolution of the gas surface density while bottom ones display that of the stellar surface density. The size of each panel is 500 pc 
by 500 pc, and its center is taken as the center of mass of the main progenitor of the final cluster. Two clusters in the center of the panel 
at 900 Myr has merged at 903 Myr. Two clusters at 903 Myr approach and merge at 906 Myr each other. The cluster in the left-bottom 
corner of the panel at 906 Myr merges with the cluster at the center by 913 Myr. 
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Fig. 5. — Time evolution of masses of hypermassive star clusters (upper panels) and their distances from the galactic center (lower 
panels). The galactic center is given by the nearest one to the cluster before the coalescence of galactic cores and by the galactic center of 
a merger remnant after that. The left and right panels show the result of runs Hpp and -ffxTi respectively. The three most massive star 
clusters are shown by the thick curves, and the others are shown by the thin curves. We show the evolution of each clusters until it sinks 
to the galactic center (the distance becomes less than 100 pc) or it is disrupted by tidal forces (Saitoh et al. 2006). Once they sink to the 
galactic center less than 100 pc or fully disrupted by the tidal force, we put circles at these points and stop tracing their evolutions. 
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Fig. 6. — Time evolution of the SFR (upper panel) and the 
bolometric luminosity (lower panel) of merging galaxies. The black 
solid and red dashed curves show results of runs Hpp and H^yt^ y 
respectively. We show the mean quantities of every 1 Myr. 



analyzed 27 samples randomly selected from the 123 I- 
band samples, which include a part of Borne's and Cui's 
samples with multiple nuclei. In their analysis, only one 
sample is classified as a ULIRG with multiple nuclei. The 
other sample includes some ULIRGs that have been al- 
ready classified as ULIRGs with multiple nuclei by I - 
band image analysis ("B orne et aLll2Q00l: ICui et al.ll2001[ ). 
The reasons for excluding them from ULIRGs with mul- 
tiple nuclei are as follows: (1) the nucleus seems to locate 
on the tidal arm and is thought to be the tidal arm or (2) 
multiple nuclei morph ology do not appear i n the iJ-band 
images. Furthermore. I Veilleux et all(|2002[ ) analyzed the 
if '-band data of 118 ULIRGs taken by the Keck observa- 
tory. Their samples also include candidates of ULIRGs 
with multiple nuclei in HST /-band data. Some of them 
are excluded since multiple cores are not apparent in K'- 
band in spite of their appearance in J-band. As a result, 
only 5 ULIRGs 4%) are classified as ULIRGs with 
multiple nuclei. 

In order to understand the difference between the re- 
sults obtained using images taken in different bands, 
we made /-, H-, and if-band images of our simulation 
data using PEGASE. We assumed that galaxies were at 
z ~ 0.1 (1" ~ 2 kpc). The point spread function was 
assumed to be Gaussian, and its dispersion was calcu- 
lated hy a — Ax/v21n2 , where Ax is the spatial res- 
olution. The spatial resolutions of /- and H-hand im- 
ages were given by the diffraction limit since we assumed 
the observations by HST. The spatial resolutions of I- 
and _ff-bands are 0".l and 0".2, respectively. For K- 
band images, the seeing determines the spatial resolu- 
tion since we assumed ground-based observations by the 
University of Hawaii 2.2 m telescope or the Keck spec- 
trosc opy. The spatial resolution, therefore, is set to be 
0".5 (H im iFaU[2002). For comparison, we also made 
the emulated images with _fC-band by James Web Space 



Telescope (JWST) or ground-based 8m telescopes with 
adaptive optics (AO). We assume that the angular reso- 
lution of these future instruments is ^ 0".05. 

Figure [8] shows /-, and K-hand images, in which 
dust extinction are not taken into account, of runs Hpp 
and 7/tt- In /-band images, multiple core-like structures 
are clearly visible, while such cores are blurred in other 
bands due to the limited resolution and the contribution 
of the luminosity of old stars. In //-band images, cores 
are connected one another and form arm-like structures. 
As a result, the cores seem to be in arms. In K-hand im- 
ages, multiple core structures are not resolved at all. As a 
result, ULIRGs seem to have only single or double cores. 
Thus, multiple cores were identified as a single or double 
cores in H- and /f'-bands. Our analysis indicates that 
ULIRGs with multiple cores can be misclassified as ones 
with single or double cores. The apparent discrepancy 
in the observational fractions of ULIRGs with multiple 
nuclei might have been caused by this difference in the 
spatial resolution. 

In the rightmost column in Fig. [SJ we show the emu- 
lated /f-band images of our simulations with the angular 
resolution of 0".05. We can clearly see the multiple bright 
sources in these /f-band images that were not visible in 
grand-base observations, since the resolution of the im- 
age is remarkably improved. We expect that, in the near 
future, observations of (U)LIRGs by JWST or ground- 
based 8m telescopes with AO will settle the problem of 
the band-to-band difference in the fraction of ULIRGs 
with multiple nuclei. 

Figure [5] shows the observed /-, //-, and /f-band im- 
ages taking dust extinction into account. Here, the cross- 
sections per H atom at H- and /T-bands are assumed to 
be aH.K = 1.0 x,10~^^ and an.K 7.0 x IQ-^^, re- 
spectively (lDrainel[2003D . We adopt — 0.5 as a fidu- 
cial value. In the /-band image, the hypermassive star 
clusters and the galactic cores are obscured because of 
dust extinction although the resolution is enough to ob- 
serve them. On the other hand, in the H- and ground- 
based /C-band images, dust extinction becomes weaker 
although the clusters and the galactic cores are not re- 
solved. In the /f-band image by JWST or ground-based 
8m telescopes with AO, the clusters and the galactic 
cores are clearly observable since not only dust extinction 
is small but also they are resolved. When dust extinc- 
tion is significant in /-band, we need the observations by 
JWST or ground-based 8m telescopes with AO. 

5. CONGLUSIONS AND DISCUSSION 

We have performed high resolution A^-body/SPH simu- 
lations of merging galaxies. We found that hypermassive 
star clusters with ^ 10* Mq form from disturbed gas 
disks in the central region kpc). In these clusters, ac- 
tive star formations take place, so that some bright core 
structures appear in the merging galaxy. The features of 
formed hypermassive star clusters agree with the obser- 
vations of ULIRGs with multiple nuclei. Although dust 
extinction may obscure the clusters after their formation, 
the clusters become optically thin within the ULIRG 
lifetime. These results indicate that one major merger 
of two spiral galaxies can explain complex structures of 
multiple nuclei in ULIRGs. Rare multiple major merger 
events on a short timescale are not necessary. 

In the high-redshift universe, gas-rich major mergers 
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Fig. 7. — Time evolution of 7-band map. The upper panels show the map without dust extinctions, and the middle panel shows the map 
taking dust extinction into account. Here, the dumpiness parameter is = 0.5. At the upper-right coner in all upper panels, simulation 
time (Myr) is shown. The bottom panels show the amount of extinctions (the unit is magnitude). 



must occur frequently. Our result suggests that these 
merging galaxies have multiple hypermassive star clus- 
ters that would be observed as bright sources in either 
near or far infrared bands. Future observations using 
e.g., JWST 01 ground-based 8m telescopes with AO are 
expected to find a number of such merging galaxies with 
multiple bright sources. 

In our simulations, we assume that main feedback 
source is Type II SNe, and radiative feedback from 
newly formed stars is not taken into account. If 
the radiative feedback is effective, gas in hypermassive 
star clusters is heated and is ejected from the clus- 
ter (|Krumholz fc Dekel l2010l) . This effect may reduce 
the cluster mass. In order to study the effect, we 
ne ed radiation-hyd r odyna mic simulations as performed 
bv lKrumholz et al.l ()2011[ ). 

Since the hypermassive star clusters are compact 
(~ IO^Mq within 20 - 30 pc), it is possible that 
star-star collisions and mergings occur in these clus- 



ters (jPortegies Zwart et al.l Il999l ). If the stellar mass 
loss in the main sequence phase is not very large, 
such merging of stars might result in the runaway 
growth of supermassive stars and t he formation of 
interniediate-mass black holes (IM BHs) (jEbisuzaki et al.l 
1200 It iPortegies Zwart et al.l [200l . These IMBHs are 
carried to the center of the galaxy together with the par- 
ent hypermassive star cluster by dynamical friction. If 
there is already a supermassive black hole at the center 
of the progenitor galaxies, the IMBH merges with the 
SMBH through dynamical friction and eccentricity evo- 
lution (Matsubayashi e t al. 2007). If there was no su- 
permassive black hole, multiple IMBHs conveyed to the 
center by the parent cluster finally merge each other to 
form the seed of a SMBH. 
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Fig. 9. — Same as Fig. |8]but taking dust extinction into account. 
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